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/  k  computer  program  has  been  put  onto  an  LSI-11  microprocessor  with  64KB  of  memory 
.  which  can  provide  accurate  ephemerides  for  GPS  (Global  Positioning  System)  satellites. 

The  satellite  dynamics  include  averaged  orbital  element  rates  due  to  J2,  tesseral 
resonances,  solar  radiation  pressure  and  third  body  perturbations  from  both  the  Moon  and 
the  Sun.  These  rates  are  first  integrated  up  to  and  across  a  satellite  pass  of  interest, 
and  a  two  point  Hermitian  interpolating  polynomial  is  established  for  each  mean  element. 
"Short  periodic  Fourier  coefficients  Hue  to  J2  and  the  Moon  and  Sun  are  next  computed,  and 
three  point  Lagrangian  interpolating  polynomials  are  established  for  them  across  the 
pass.  .Both  sets  of  interpolating  polynomials  are  finally  used  to  provide  osculating 
^orbital  elements  at  arbitrary  times  during  the  pass. 

The  modular  form  of  the~computer  program  and_'the  sequential  structure  of  the  theory 
mean  that  no  large  portion  of  the  program  needs  to  exist  in  core  at  once.  The  use  of 
overlays,  then,  allows  the  small  computer  to  handle  the  large  program;  - 

The  computation  of  the  interpolating  coefficients  for  orbital  element  prediction 
during  a  satellite  pass  seven  days  beyond  the  epoch  time  takes  about  100  seconds  in  real 
time.  Once  the  interpolating  coefficients  are  computed  the  state  vector  of  the  satellite 
at  some  desired  time  during  the  satellite  pass  can  be  obtained  in  less  than  1  second  of 
^  real  time.  Comparisons  with  comparable  runs  with  the  Draper  Lab  GTDS  R&D  program  have 
shown  discrepancies  in  position  less  than  20  meters. 

—  "t' ^^JhpS  ^computer  program  includes  an  analytical  Lunar/Solar  eohemeris  so  it  is  self- 
contained  except  for  input  mean  orbital  elements.  Partial  derivatives  have,  been  im¬ 
plemented  which  will  give  the  capability  to  fit  observations  of  the  satellites  and  to 
consequently  obtain  the  necessary  mean  elements. 

The  program  can  be  modified  quite  easily  to  handle  synchronous  satellites  by  modify¬ 
ing  the  subroutine  modules  for  tesseral  resonant  perturbations  and  lunar-solar  short- 
. periodics.  With  the  present  overlay  scheme,  considerable  expansion  of  the  program  is 
possible  to  obtain  more  accuracy  and  versatility 

♦Research  at  MIT  supported  by  U.S.  Air  Force  Geophysics  Laboratory,  Geodesy  and  Gravity 
Branch,  under  contract  F19628-82-K-0002 .  \ 
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A  computer  program  has  been  put  onto  an  LSI-11  micro¬ 
processor  with  64KB  of  memory  which  can  provide  accurate 
ephemerides  for  GPS  (Global  Positioning  System)  satellites. 

The  satellite  dynamics  include  averaged  orbital  element 
rates  due  to  ,  tesseral  resonances,  solar  radiation  pressure 
and  third  body  perturbations  from  both  the  Moon  and  the  Sun. 

These  rates  are  first  integrated  up  to  and  across  a  satellite 
pass  of  interest,  and  a  two  point  Hernitian  interpolating 
polynomial  is  established  for  each  mean  element.  Short  periodic 
Fourier  coefficients  due  to  2nd  the  Moon  and  Sun  are  next 
computed,  and  three  point  Lagrar.gian  interpolating  polynomials 
are  established  for  them  across  the  pass.  Both  sets  of  inter¬ 
polating  polynomials  are  finally  used  to  provide  osculating 
orbital  elements  at  arbitrary  times  during  the  pass. 
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The  modular  form  of  the  computer  program  and  the  sequential 
structure  of  the  theory  mean  that  no  large  portion  of  the  program 
needs  to  reside  in  core  at  once.  The  use  of  overlays,  then, 
allows  the  small  computer  to  handle  the  large  program. 

The  computation  of  the  interpolating  coefficients  for 
orbital  element  prediction  curing  a  satellite  pass  seven  days 
beyond  the  epoch  time  tahes  about  100  seconds  in  real  time.  Once 
the  interpolating  coefficients  are  computed  the  state  vector  of 
the  satellite  at  some  desired  tire  during  the  satellite  pass  can 
be  obtained  in  less  than  1  second  of  real  time.  Comparisons  with 
comparable  runs  with  the  Draper  Lab  GTDS  R&D  program  have  shewn 
discrepancies  in  position  less  than  20  meters. 

The  computer  program  includes  an  analytical  Lunar/Solar 
ephem.eris  so  it  is  sel f-cor.ta : r.ec  except  for  inout  mean  orbital 
elements.  Partial  derivatives  have  been  implemented  which  will 
give  the  capability  to  fit  observations  of  the  satellites  and  to 
consequently  obtain  the  necessary  mean  elements. 

The  program  can  be  codified  quite  easily  to  handle  syn¬ 
chronous  satellites  by  modifying  the  subroutine  modules  for 
tesseral  resonant  perturbations  and  lur.ar-solar  shor  t-per  iodics . 
With  the  present  overlay  scheme,  considerable  expansion  of  the 
program  is  possible  to  obtain  mors  accuracy  and  versatility. 


INTRODUCTION 


The  motivation  for  this  work  came  as  a  result  of  desiring 
accurate  ephemerides  for  GPS  (Global  Positioning  System) 
Satellites  which  could  be  computed  on  an  LSI-11  microprocessor 
with  64  KB  of  memory.  The  satellite  ephemerides  are  required  for 
MITES  data  reduction  which  is  done  on  the  LSI  computer.  MITES  is 
a  system  of  miniature  radio  interferometer  terminals  for  the 
measurement  of  baseline  vectors  on  the  ground  by  means  of 
interferometric  observations  of  C-PS  satellites  (Ref.  1). 

A  model  was  needed  which  would  run  reasonably  fast  and  of 
course  would  fit  into  the  small  computer.  Although  not  attained 
yet,  our  eventual  accuracy  goals  are  at  the  one  part  per  million 
level  or  better. 

The  needs  of  the  MITES  project  seemed  to  best  be  satisfied 
with  the  GTDS  semianalyt ical  orbit  model  used  on  the  AMDAHL 
computer  at  Draper  Labs.  The  advantage  of  using  this  type  of 
model  over  present  analytical  (General  Perturbations)  theories  is 
the  generality  of  the  force  models  which  can  be  incorporated  and 
the  ease  with  which  these  models  may  be  included  or  replaced. 

The  latter  is  important  for  tailoring  the  full  scale  GTDS  version 
to  the  LSI. 

The  tailoring  study  will  be  presented  in  mere  detail 
following  a  general  description  of  the  model.  Following  that  the 
various  force  models  will  be  referenced  and  discussed.  Then  a 
discussion  of  the  computer  program  is  given  followed  by  the 
results  of  comparison  with  mere  complete  orbit  theories. 


GENERAL  DESCRIPTION  OF  THE  MODEL 

V7e  begin  by  giving  a  general  overview  of  the  rcdel. 

Instead  of  the  more  traditional  Keplerian  elements  the  model 
uses  equinoctial  elements  which  are  expressed  in  terms  of  the 
Keplerian  elements  viz 
a  =  a 

h  =  e  sin(d)  +  fl) 
k  =  e  cos(u  +  n) 
p  =  tan ( i/2 ) sin (a) 
q  =  tan ( i/2) cos (a) 
x  =  m  +  b)  +  a 

where  a  =  semimajor  axis,  e  =  eccentricity,  i  =  inclination,  M  = 
mean  anomaly,  u  =  argument  of  pericenter,  and  a  =  longitude  of 
ascending  node. 

These  elements  vary  secularly  with  short  periodic  effects 
superimposed.  The  osculating  value  of  an  element  at  some  time  is 
therefore  the  sum  of  the  mean  value  of  the  element  and  the  short 
periodic  contributions.  One  method  of  computing  the  osculating 
values  is  with  a  three  part  model  which  includes: 

1)  the  computation  of  the  mean  values  of  the  elements; 

2)  the  computation  of  the  coefficients  fer  the  Fourier 
series  expansions  of  the  short  periodics  of  the  elements;  and 

3)  the  combination  of  the  results  of  the  two . 

This  is  accomplished  as  follows: 

The  mean  values  are  determined  from  averaged  rate  equations 
which  are  integrated  with  a  fourth  order  Pur.ge  Kutta  technique. 
The  rate  equation  for  each  element  is  the  s<jm  of  individual  rate 
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equations  due  to  the  various  perturbations.  The  rate  equations 
are  integrated  from  a  specific  epoch  (when  initial  mean  elements 
are  given)  up  to  the  beginning  of  a  satellite  pass.  The  typical 
step  size  of  the  Runge  Kutta  integrations  is  1/2  cay.  The  rate 
equations  are  then  integrated  in  one  step  to  the  end  of  the  pass. 
Using  the  resultant  equinoctial  elements  at  the  beginning  and  end 
of  the  pass,  a  two  point  Hermitian  interpolating  polynomial  is 
determined  for  each  element.  This  polynomial  provides  the 
averaged  value  of  the  element  at  any  time  during  the  pass. 

The  short  periodic  Fourier  coefficients  are  computed  at 
three  epochs  across  the  pass  -  the  beginning,  middle  and  end. 

The  Hermitian  interpolating  polynomials  are  used  to  obtain  the 
mean  elements  required  for  the  evaluation  of  the  short  periodic 
coefficients.  A  sine  and  cosine  coefficient  for  each  of  four 
frequencies  of  X  is  computed  as  a  sum  over  all  contributing 
perturbations.  A  three-point  Lagrangian  interpolating  polynomial 
is  then  determined  for  each  coefficient. 

Finally,  there  is  the  construction  of  the  osculating  equi¬ 
noctial  elements  at  a  desired  time  during  the  pass.  This  is  done 
simply  by: 

a)  obtaining  the  mean  elements  from  the  Hermitian 
interpolators; 

b)  obtaining  the  short  periodic  coefficients  from  the 
Lagrangian  interpolators; 

c)  evaluating  the  Fourier  series  with  the  short  periodic 
coefficients  and  the  mean  values  of  X  to  obtain  the  short 
periodic  perturbations  in  the  elements;  and 


d)  adding  the  mean  elements  to  the  short  periodic  variations 
to  obtain  the  osculatina  values. 


TAILORING  OF  THE  MODEL  FOR  A  SMALL  COMPUTER 

The  complete  GTDS  semianalytical  theory  provides  10  cm  to  1 
meter  accuracy  in  comparison  with  a  precision  (Cowell)  numerical 
integration  of  the  GPS  orbit  (Ref.  2).  The  precision  integration 
was  performed  using  a  GTDS  program  as  modified  at  Draper  (Ref.  3) 
and  included  all  the  terms  significant  for  GPS: 

-  8  x  8  gravity  field 

-  lunar-solar  point  masses 

-  solar  radiation  pressure. 

Designing  a  tailored  theory  first  means  specifying 
truncations  for  the  rate  equations  and  for  the  short-periodic 
models.  The  design  of  the  tailored  theory  is  greatly  aided  by 
output  reports,  from  the  GTDS  semianalytical  program,  of  the  mean 
element  rates  and  the  short  periodic  coefficients.  Summaries  of 
these  reports  are  given  in  Tables  I  (for  the  mean  element  rates) 
and  II  (for  the  short-periodic  coefficients). 

In  analyzing  the  data  of  Table  I,  the  following  error  bounds 
were  adopted: 

(At)  (a5)  <  20  meters 
(a)(At)(iH)  <  20  meters 
(a)  (tt)  (Ale)  <  20  meters 
(a) (At) (ip)  <  20  meters 
(a)  (At)  (Aa)  <  20  meters 
(a) (At) (aI)  <  20  meters 


Terra 

Maanit 

J2 

2  km 

Lunar-Solar 

200  meters 

Point  Masses 

Tesseral  M-Dailies 

20  meters 

Tesseral  Linear 

15  meters 

Combination  Terras 

Solar  Radiation 

15  meters 

Pressure 

J2 

1 5  cm 

where  a,  the  semimajor  axis  of  the  GPS  orbit  is  in  meters  and  it 


is  an  approximation  interval  in  seconds.  These  bounds  relate,  in 
a  very  rudimentary  wav,  to  bounds  in  the  radial  (h  and  k) ,  cross 
track  (p  and  q)  and  along  track  ( X )  errors.  Assuming  a  =  25559 
km  and  At  =  14  days  gives  an  element  rate  error  bound 

Aa.  <  6  x  10~13. 

1  • 

Based  on  this  bound,  we  have  currently  adopted  all  of  the 

averaged  equation  of  motion  models  in  Table  I  except  for  those 

2 

for  Jj.  Higher  degree  rcnals  (J^,  ,  etc.)  were  eliminated 

either  due  to  the  near  critical  inclination,  the  high  altitude 
and  small  Rg/a  (where  R^  is  the  Earth  radius)  or  the  small  eccen¬ 
tricity  of  the  GPS  orbit.  For  the  tesseral  resonance,  it  was 
known  from  a  previous  study  (Pef.  4)  that  rhe  (C,S)3  2  terms  were 
dominant.  Therefore,  the  nominal  tesseral  resonance  model  should 
include  at  least  the  (C,S)3  2  terms.  Reference  4  and  further 
study,  which  will  be  discussed  below,  also  indicates  that  the 
(C,S)2  2  an<^  (C,S)^  ^  terms  have  to  be  included.  Full  recursive 
models  for  the  lunar-solar  point  mass  terms  (Ref.  5)  and  a 
numerical  averaging  algorithm  for  solar  radiation  pressure  (Ref. 
6)  were  adopted  for  the  testing. 

Based  on  Table  II  and  the  20  meter  requirement,  it  was 
decided  to  retain  the  J2  zonal  short  periodic®,  the  ?2  third-body 
short  periodics  and  the  tesseral  m-dailies.  Further,  it  was 
decided  not  to  include  any  vreak-Time-Oepenbent  (WT3)  (Ref.  7) 
corrections  to  the  third-body  short-periods.  The  WTD  terms  are 
budgeted  at  10  meters  for  the  GPS  case.  it  is  estimated  that  the 
Pj  third-body  short  periodics  are  15  meters  for  this  case. 
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To  examine  the  assumptions  of  the  tailored  theory,  tests 
were  made  against  the  GTDS  precision  numerical  integration.  The 
procedure  is  as  follows: 

1.  generate  a  truth  ephemeris  with  the  precise  nu-erical 
integration  (Cowell)  with  its  complete  force  models; 

2.  least  squares  fit  the  tailored. semianalytical  satellite 
thoery  to  the  truth  ephemeris; 

3.  generate  statistics  and  plots  for  the  residuals  between 
the  semianalytical  ephemeris  (using  the  adjusted  intitial 
conditions  from  the  least  squares  fit)  and  the  truth  echemeris 
constructed  in  step  1. 

The  major  results  are  given  in  Table  III.  The  m-daily 
tesseral  short  periodica  were  not  included  in  this  testing. 

AVERAGED  RATES  AND  SHORT  PERICDICS 

With  the  exception  of  solar  radiation  pressure,  the  force 
models  and  development  of  the  averaged  rate  equations  and  short 
periodic  expansions  for  the  perturbations  have  been  documented  in 
the  literature.  With  regard  to  these  then,  we  will  onlv  show  a 
table  of  references  (Table  IV)  and  make  a  few  specific  comments 
about  their  implementation  cn  the  LSI.  A  simplified  analytical 
model  has  been  implemented  for  the  solar  radiation  pressure 
averaged  rate  equations.  As  this  contribution  to  the  theory  has 
not  been  documented  in  the  literature,  the  -ath.e-atical 
development  will  be  briefly  summarized. 

Two  forms  of  the  J2  averaged  rate  equations  were  constructed 
for  the  LSI  model  -  a  truncated  form  and  a  closed  form  in  the 
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eccentricity.  The  forms  afe  nearly  identical  except  that  the 
closed  form  has  an  additional  factor  of  X  to  a  power  where  X  * 

1/ .  Little  is  gained  with  regard  to  accuracy  by  using 
the  closed  form,  but  on  the  other  hand  little  time  is  lost 
computationally.  Therefore  the  closed  form  has  been  adopted. 

Initially  it  was  felt  that  due  to  the  low  eccentricity  of 
the  GPS  satellites  (e  <  0.009}  the  Fourier  expansions  of  the 
short  periodic  generator  need  include  just  3  frequencies  of  I. 
Early  testing  of  the  LSI  model  against  the  GTDS  model  for  a  GPS 
satellite  with  e  3  0.003  seemed  to  indicate  that  this  was  a  good 
assumption.  In  tests  involving  a  GPS  satellite  with  e  *  0.009, 
however,  discrepancies  larger  than  expected  occurred.  These  were 
eventually  eliminated  by  including  frequency  terms  in  the 
Fourier  expansions. 

The  rate  equations  that  were  constructed  for  the  tesseral 
resonances  were  all  truncated  to  0th  order  in  the  eccentricity. 

The  tesseral  ir.-daily  model  has  not  been  fully  implemented  in 
the  LSI  yet  and  so  the  results  which  will  be  presented  below  do 
not  contain  their  effect. 

The  computation  of  the  averaged  differential  equations  and 
short  periodics  for  the  lunar-solar  perturbations  require 
ephemerides  for  the  Moon  and  Sun.  For  the  LSI  model  we  decided 
to  utilize  the  low  precision  Fourier  series  formulae  of  T.  C.  van 
Flandern  and  K.  F.  Pulkkinen  of  the  DSh’O  for  the  geocentric 
positions  of  the  Moon  and  Sun  (?ef.  12}.  The  original  code  is 
very  lengthy,  requiring  >  24F3  for  the  Moon  and  »  4.3K3  the  Sun. 
Therefore,  for  the  initial  implementation  of  the  OSSO  lunar/solar 
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ephemeris  on  the  LSI,  we  used  a  truncated  form  of  the  formulae 
due  to  Belastock  (Ref.  13)  which  deleted  terms  of  order  x  10-<* 
and  x  10-5.  In  a  study  similar  to  that  which  was  used  to  prepare 
Table  III,  Belastock  showed  an  error  of  about  20  meters  due  to 
this  truncation.  Eventually  we  were  able  to  fit  the  non- 
truncated  series  expressions  into  the  LSI  (by  breaking  the 
subroutine  for  the  Moon  into  four  sections  and  overlaying  each 
part  against  the  other),  and  eliminated  this  source  of  error. 

The  following  briefly  reviews  the  development  of  the 
simplified  analytical  model  for  solar  radiation  pressure.  Ke 
should  point  out  that  the  only  new  aspect  of  this  model  is  its 
development  in  terms  of  equinoctial  elements. 

For  special  perturbation  integrators  such  as  the  G7DS 
system,  the  solar  radiation  pressure  is  modeled  with  the  formula 
(Ref.  14) . 

2  CRA 
=  vPeR2 

SR  S  sun  _  ,3 

®  I Kyg 1 

where 


eclipse  factor, 

{mean  solar  flux  at  1  AE'/fspeed  of  light), 

1  AO, 

1  surface  reflectivity, 
effective  area. 


mass  , 

r  -  Pg  where 

r  *  position  vector  of  the  spacecraft  in  inertial 


coordinates 


7  V  V /  / 
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Rg  =  position  vector  of  the  Sun  in  inertial 


coord inates . 


In  the  following  the  satellite  is  assumed  to  always  be  in 


sunlight  so  that  v  =  1  and  the  factor  P_R  -  is  constant 

S  sun  m 


Wi  th 


2 

T  =  vP  R  - 

S  sun  m 


then 


lBvsl 


The  geometrical  situation  is  similar  to  a  third  bodv  mass 


perturbation  where  the  Earth-center  to  satellite  distance  is  very 


small  relative  to  the  Earth-center  to  third  bodv  distance.  This 


leads  us  to  write  the  equivalent  disturbing  potential  (Ref.  15) 


where  A  as  usual  is  written 


A  =  Rs[l  +  <£->2-  jr  cost]  1/2 
S  S 


and  *  is  the  angle  between  r  and  p.g.  with  the  Legendre  ex?ansi< 


(Ref.  16),  we  obtain  the  disturbing  potential 


11  °SR  *  57  r,  (S7>  pn(cos*' 

S  n=  l  s 


The  limitations  of  this  equation  are: 


1.  the  surface  reflectivity  is  constant; 


2.  the  cross-sectional  area  perpend icular  to  the  flux  i: 


constant; 


3.  the  satellite  is  not  eclipsed. 


*  '-Jk  "-I 


The  third  assumption  means  that  during  eclipsing  periods  the 
potential  cannot  strictly  account  for  solar  radiation  pressure. 
For  GPS  satellites  those  periods  occur  for  4-5  weeks  twice  a 


year 


The  above  assumptions  and  the  smallness  of  r/R  succests 


that  equation  1)  be  written: 


2)  U, 


-Trcos*  -T(Rs‘r) 


Four  more  stages  of  algebra  lead  to  the  desired  differential 
equations : 

1.  substitution  of  Rg  =  of  +  6g  +  rw  and  r  =  (rcosL)f  + 
(rsinL)g  into  equation  2)  where  a,  B  and  y  are  the  direction 
cosines  of  the  Sun  relative  to  the  equinoctial  f,  g,  w  reference 
system; 

2.  averaging  over  the  mean  longitude  L; 

3.  taking  partial  derivatives  of  Ug.  with  respect  to  a,  h, 

k,  a  and  3; 

4.  substituting  these  resultant  expressions  into  ecuation 
(3-71)  of  P.ef .  17). 

The  final  expressions  are: 


=  0 


3T  =  Va*t®  +  ]S(hq  '  K?)- 


=  -Va[|s  +  £g(hq  -  <?>  I 


a?  =  va^h^ 


at  *  Va(ix5,kT 

at  *  var-|(«k  ♦  eh)  +  -  ksi  ♦  i_(hq  -  kPn 


where 


A  =  na  ,  n  =  mean  motion 

B  =  A  -  h2  -  k2 

C  =  1  +  p2  +  q2. 

PARTIAL  DERIVATIVES 

The  purpose  of  computing  partial  derivatives  in  the  model  is 
to  determine  from  data,  in  a  least  souares  sense,  the  mean 
elements  at  epoch.  The  data  can  be  in  two  forms: 

-  observations 

-  satellite  state  vectors  generated  by  numerical  integra¬ 
tion. 

If  one  if  fitting  to  numerically  integrated  data  the  general 
expression  for  the  partial  is 

9x,  x  3 x ,  X  .  3fosc.  elm.)  f  .  3 (mean  elm.)f 

3  (mean  elm.)”  1 1  ~  3  (osc.  e.m. ,  I  a  (mean  elm.  j  1 1  3  (mean  iim.) 

o  i 

The  first  matrix  on  the  right  of  this  equation  is  derived  from 
two  body  equations  (Ref.  7),  the  second  matrix  for  our  purposes 
is  the  identity  matrix  and  the  third  matrix  has  variational 
equations  (Ref.  7)  which  must  be  integrated  via  Runge  Kutta  at 
the  same  time  that  the  averaged  rete  equations  are  integrated  in 
the  model. 

If  one  is  fitting  to  observations  one  has 

3  fobs.)  a  »  fobs.)  .  3  x  ,  x 

3 (mean  elm.)..  *  3(mean  elm.).  * 


THE  COMPUTER  PROGRAM 

The  LSI-11  has  64  KB  of  main  memory.  Approximately  8  KB  of 
this  is  occupied  by  system  routines  leaving  56  KB  for  the 
program  to  be  executed.  Fortunately,  due  to  the  overlay 
capability  of  the  LSI,  this  is  not  the  upper  limit  to  the  size  of 
a  program.  A  program  using  overlays  can  be  considerably  larger 
than  would  normally  fit  into  the  central  memory  since  portions  of 
the  program  reside  on  disk. 

The  basic  requirement  for  defining  an  overlay  structure  is 
that  only  relatively  small  parts  of  the  program  be  required  for 
execution  at  any  one  time.  An  attractive  feature  of  the  semi- 
analytical  model  for  implementation  on  a  small  computer  is  the 
sequential  nature  of  its  computations  which  makes  it  quite 
suitable  for  overlaying.  Figures  1,  2,  3  and  4  show  the 
subroutine  structure  of  the  program.  Each  diagram  is  a  complete 
path  of  computations  requiring  only  the  end  results  of  the 
previous  diagrams.  Kithin  each  diagram,  as  well,  a  distinct 
sequential  structure  is  clear. 

On  the  LSI  the  overlay  structure  which  we  set  up  consists  of 

a  root  and  three  overlay  regions.  The  root  contains  those  parts 

of  the  program  that  must  always  reside  in  memory  which  in  our 

case  is  the  main  program  and  a  few  global  subroutines.  Each  of 

the  three  overlay  regions  contain  a  number  of  seg-ents.  Each 

segment  consists  of  a  group  of  subroutines  which  are  independent 

• 

of  subroutines  in  other  segments  of  that  region,  i.e.  from  one 
segment  of  the  region  one  can’t  call  a  routine  in  another  segment 
of  the  same  region.  This  is  important  since  only  one  segment 
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Fig.  2  Continued 


from  a  region  can  be  resident  in  central  memorv  at  a  time. 
Otherwise  the  return  path  to  the  main  orocram  would  be  destroyed. 
It  is  necessary  though  to  make  use  of  the  fact  that  segments  in 
one  region  can  make  calls  to  routines  in  segments  of  other 
regions.  This  is  why  three  regions  were  established.  Figure  5 
shows  how  the  overlay  regions  are  constructed  and  interlinked. 

In  going  to  a  small  computer  one  expects  a  longer  execution 
time  for  the  program.  The  timing  of  the  program's  execution  is 
separated  into  two  parts:  1)  the  time  required  to  set  up  all  of 
the  interpolation  coefficients  -  =10C  seconds;  and  2)  the  time 
required  to  compute  a  satellite  position  from  them  -  =1  second. 

RESULTS 

During  the  development  cf  the  LSI  model  two  types  of 
comparisons  were  mace  for  debugging  the  program  end  for  checking 

the  error  level 

-  comparison  with  the  GTDS  semianalytical  model 

-  comparison  with  numerical  integration. 

Most  of  the  comparisons  with  the  GTD3  model  were  for 
debugging  purposes  since  the  averaged  rate  equation  and  short 
periodic  computations  could  be  individually  checked  for  each 
force  model.  The  most  significant  comparison  was  made  against  a 
7  day  GTDS  run  where  both  models  used  the  sa-e  initial  mean 
elements.  Agreement  in  position  during  the  7  day  span  was  on  the 
order  of  20  meters. 

For  examining  the  error  level  of  the  LSI  model,  comparisons 


were  mostly  made  against  numerical  integration 


Two  oossible 


numerical  integration  programs  existed  for  the  purpose  of 
comparison  -  the  GTDS  system  already  mentioned  and  PEP  (Pef. 
which  is  used  at  MIT  and  is,  as  GTDS,  a  full  blown  model.  As 
transfer  of  comparison  files  from  the  PEP  computer  program  to 
LSI  was  more  convenient  than  from  the  AMDAHL,  PEP  was  used. 

Given  a  set  of  high  precision  orbital  elements  at  epoch. 


18) 

the 

a 


Sv 


corresponding  set  of  initial  mean  elements  must  be  determined  for 
the  semianalytical  model.  It  is  important  to  obtain  these  mean 
elements  accurately  since  errors  in  them  lead  to  a  secularly 
increasing  prediction  error.  A  least  squares  differential 
correction  algorithm  is  used  that  solves  for  the  eooch  mean 
elements  based  on  the  fit  of  a  semianalytical  trajectory  to  the 


ephemeris  output  of  a  high  crecisicn  Cowell  integration.  The  a 


priori  estimates  of  the  mean  elements  at  epoch  are  obtained  from 
a  single  point  conversion  method.  Here  the  short  oeriouics  of 
the  model  are  computed  at  the  epoch  and  subtracted  from  the 
osculating  elements.  PEP  was  used  to  generate  the  high  orecision 
ephemeris  for  a  fit  span  of  36  hours  with  a  22.5  minute  spacing. 

Osina  the  mean  elements  thus  obtained,  a  semianalytical 
trajectory  was  generated  again  -  this  time  for  ccmoarison  pur¬ 
poses  with  the  PEP  ephemeris.  Fias.  6,  7  ar.c  3  show  respectively 
radial,  along  track  and  cross  track  cc-oariscns  with  PEP  over  the 
36  hour  fit  span.  Figs.  9,  10  and  11  show  comoarisons  over  the 
36  hour  fit  span  and  the  36  hours  beyond  the  fit  scan.  The 
latter  three  figures  show  an  expanding  envelope  of  the  residuals 
with  tine.  This  is  due  to  various  model  differences.  It  appears 


that  the  error  level  of  the  LSI  eodel  is  on  the  order  of  2  pom 
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Fig.  10  LSI  vs  PEP  in  along  track,  the  first  36 
hours  are  over  the  fit  span  and  the  last 
36  hours  are  over  the  predict  span. 
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Fig.  11  LSI  vs  PEP  in  cross  track,  the  first  36 
hours  are  over  the  fit  span  and  the  last 
36  hours  are  over  the  predict  span. 


FUTURE  WORK 


With  the  present  overlay  structure  we  have  the  flexibility 
of  improving  the  accuracy  of  the  model  by  incorporating  some  of 
the  force  models  we  have  so  far  omitted.  We  would  eventually 
like  to  handle  solar  radiation  pressure  .in  a  more  complete  manner 
since  this  is  an  important  effect  for  the  high  altitude  GPS 
satellites.  The  program  can  be  quite  easily  suited  to  syn¬ 
chronous  satellites  by  modifying  the  subroutine  modules  for 
tesseral  resonant  perturbations  and  the  lunar-solar  short- 
periodics.  Finally,  we  have  immediate  plans  to  fit  real 
observations  with  the  LSI  model  in  order  to  determine  initial 
mean  elements. 
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